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Abstract 



in gas filled hollow-core photonic crystal fibers that includes a plasma induced nonlinearity. As 
anticipated for a simpler model and using a perturbation analysis, there are indeed stationary soli- 
tary waves that accelerate and self-shift to higher frequencies. However, if the plasma nonlinearity 



> 

o , 

f^ I strength is large or the pulse amplitudes are small, the solutions have distinguished long tails and 



decay as they propagate. 
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I. INTRODUCTION 



Hollow-core photonic crystal fibers (HC-PCFs) can exhibit very interesting properties, 
such as relatively low loss, low group velocity dispersion and high confinement of light in 
the core [l|, 1^, while also allowing new nonlinear phenomena associated with the interaction 
of light and matter filling these fibers. Lately, these HC-PCFs have been filled with gases 
for purposes of enhancing Raman scattering if a Raman-active gas is used 3|], or further 
controlling the total dispersion of the fiber by varying the gas pressure j^. Furthermore, it 
has also been shown that few juJ or even picoJ 5|, Isl energy optical pulses are sufficient to 
ionize the gas and produce a plasma, leading to new nonlinear effects such as the blueshifting 
of the central wavelength of the pulses [a, [3]. Despite the fact that the soliton shift to higher 
xequencies has also been reported in other contexts, such as in a line-defect waveguide 
8| and in tapered solid-core photonic crystal fibers [9j, the existence of a blueshift in a 
Raman-active gas opens new exciting opportunities of controlling the soliton dynamics by 
two competing processes, one leading to a redshift, usually known as soliton self-frequency 
shift (SSFS) caused by intrapulse Raman scattering (IRS) |lO|, and the other to a blueshift. 

Traditionally, the interaction between light and matter has been studied using compu- 
tationally demanding methods based on models for the full electric field of the pulse 11 1 
but, recently, Saleh et al. presented a model that describes pulse propagation in hollow-core 
photonic crystal fibers filled with a gas as a pair of coupled equations for the electric field 
envelope and ionization fraction 12|]. This model, which neglects losses and results from a 
linearization of the tunneling model for pulse intensities close to the threshold intensity, has 
proved to be amenable to the application of both numerical and analytical techniques. In 
effect, by using a perturbation approach the occurrence of the blueshift effect has already 
been adequately predicted |12l |. 

In this work, we present a thorough study of accelerating solitons in gas-filled HC-PCF, 
extending the results in 12|. We start with the model proposed by Saleh et al. |12|, use 
an accelerating self-similarity variable to obtain an ordinary differential equation (ODE) 
to which we apply a perturbation approach and solve using a shooting procedure. In this 
analysis, we have considered the exact solution for the ionization term and our results apply 
to both zero and nonzero threshold intensities. The dependence on the model parameters, 
namely, the plasma and Raman strengths, the intensity threshold and the pulse peak value 



are studied in detail. 



II. SELF-SIMILARITY VARIABLE AND PERTURBATION APPROACH 



As mentioned in the introduction, here we will follow Saleh et al. |12| and start with the 
following coupled equations 

-^ = {a/A,s){nT - n,)A\^fe{A\^|Jf) (2) 

where ip{z,t) is the optical field envelope in units of square root of power, z is the distance 
along the fiber, t is the time in a reference frame moving with the group velocity at the central 
frequency Uq, P2 is the group velocity dispersion parameter, 7 is the nonlinear parameter, 
tji is the Raman parameter, c is the vacuum speed of light, ko = ojq/c, Up = [e^ne/(eome)]^/^ 
is the plasma frequency associated with an electron density ne{t), e and mf. are the electron 
charge and mass, respectively, eo is the vacuum permitivity and Aes is the effective mode 
area. The plasma-induced nonlinearity only occurs for intensities above the threshold inten- 
sity Jth = iV'lthMeflf, SO that A|?/'p = |'?/'p — It/'I^j^ and 9 is the Heaviside step function, ut is 
the total number density of ionizable atoms, associated with the maximum plasma frequency 
cot = (e^riT / {eoiTT'e)Y^'^ and a is the photoionization cross-section. This model assumes that 
the recombination time is longer than the pulse and neglects the ionization induced loss that 
is small especially for pulses whose maximum is barely above the threshold. The equation 
([2]) may be solved exactly and after an adimensionalization we obtain 



I + ^0 + krg - rn{\q\\<l - <Pt, (1 - e--/^^^^l^l^«(^l^l^)-') = (3) 



where q = {'jzoY^'^^, ^ = z/zq, t = t/to, tr = tn/to, (px = \koZQ{ujT / (^of , o = crto/{A^s'yzo), 
where zq = to/l/^al is the called the dispersion length and to is an arbitrary time chosen 
similar to the pulse duration. Motivated by the observation of blueshifting of the pulses and 
previous perturbation approaches, we used an accelerating variable T = t + |^^ + bC, and 
solutions of the form g(^, r) = F(T) exp(i6'(^, T)), with F and 9 real, and obtained an ODE 
for F 



F" + aTF -DF + 2F^ - AtrF'F' - 20tF (l - e-'^^-^o AF^e(AF2)dT' j ^ q^ 

and the following expression for the phase 

0{^, T) = - (^e + 6) T + i(D + b')^ + ^bae + ^a'e + E (4) 

where D and E are arbitrary constants. In order to reduce the number of parameters, we 
introduced the following change of variables -P(C) = <^F{T) and T = aC,., with which the 
ODE for P(C) reads 

P" + aCP -CP + 2P3 _ ^^p2p/ _ ^^p (^^ _ g- /i^ Ap2e(Ap2)dc' j = (5) 

with a = aa^, 7ij = Arn/a and 7p = 20^0"^ • If we further define 

7p = X, 7i? = ^X 

where R = 'Jr/'^p, it will permit the application of a perturbation approach simultaneously 
to the two terms, namely, Raman and plasma, as long as x is small. 

Hence, we have used a perturbation approach around the ODE associated with the nonlin- 
ear Schrodinger equation (NLSE) whose results are valuable by themselves if the additional 
terms are small, but that also serve as first estimates for our shooting method. Hence, we 
consider expansions for P and a in powers of x such that 

p = GiC-Co) + xPi{C) + ---, 

where G(C - Co) = VCsech{VC{C - Co)) and 

a = x"! + ■ ■ ■ 
and introduce them in ([5]). To first order, we obtain 

P{' - CPi + 6G^Pi = -aiCG + RG^C 

_ G (1 - e-^--(^'-^'h)e(G2-nl)'iC'A 

The left member of the last equation is obeyed by G', so that, the solvability condition is 

(GG'dt = R G'^{G'fdC 
-00 J —00 

+ r GG' (1 - e-/-o.(c^^-^.i)e(G^-^4)'^C'^A ^ 
J —00 

which gives 
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+ J^fl_e-2v/^^^+2P.lCi) (6) 

where C' = C ~ Co and Ci is the instant at which (^(±(^1) = Pth- The integral in the 
last expression may be written in closed form as a series but here we solve it numerically. 
Nevertheless, for Pth = 0, the integral is easily solved analytically such that a\ simplifies to 
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In the limit of small C, this equation reduces to ai = —^RC^ + |C Moreover, a graphical 
inspection of expression ([7]) shows that it is always below the curve |C and in fact it tends 
to 1 as C increases. On the other hand, whenever the peak intensity is close to the intensity 
threshold, i.e., y/C ~ Pth, we may expand the exponentials in equation (|6]) up to first order 
and obtain 

a, = -^RC + -^{C - Plfl\ (8) 

15 sVC 

Note that both this expression and the approximate expression for the acceleration for 

small C and Pth = exhibit C^ and C dependencies which are associated, respectively, to 

Raman and plasma effects. Such dependencies imply that the plasma effect is expected to 

dominate for small peak amplitude pulses, with the acceleration taking positive values which 

are proportional to the square of the peak amplitude. Conversely, as the peak amplitude 

increases, the soliton trajectory should be mainly controlled by Raman effect, which leads 

to a negative acceleration that is dependent on the forth power of the peak amplitude. 



III. PULSE PROFILES AND ACCELERATIONS 

We then used a shooting method to obtain the pulse solutions of equation (E]) and respec- 
tive accelerations. For this purpose, we first analyse the asymptotical form of this equation 
for pulse solutions that vanish at the limits C, — > ±00 which is given by 

P" + K -C-xAoo)P = 0, 



where 

A = / ° if C ^ -oo 

that, using z = —a^^^C + a^^/^(C + x^oo), niay be transformed to an Airy equation: 

p" _ ^p = 0. 

This resuh anticipates that the pulse solutions with a > have tails that are exponentially 
decreasing as C — ?■ — oo as the Airy solution Ai(z) for z — t- oo and may have tails that are also 
exponentially decreasing as ( increases as the solution Bi{z) for z ^ 0^ but that eventually 
exhibit Airy oscillations (even if very, very small). For a < 0, the contrary is true. 

Considering those asymptotic behavior, we designed our shooting procedure as follows. 
First, we have fixed the acceleration a and, starting from the left tail, we integrated forward 
using initial conditions that conform with the corresponding Ai(2;). The actual location of 
the pulse in the ( axis may be estimated using the perturbation approach described in the 
previous section but, since it is not very far from C = 0, the first estimate for the left tail 
location Cminus may be obtained as if -P(C) ~ G{(). Then, the shooting procedure checks if P 
and P' are already very small at some point in the right tail, and improve the starting ( in 
order to obtain the actual pulse profile for the chosen acceleration. Therefore, this procedure 
allow us to establish the relationship between the acceleration and the pulse characteristics, 
namely, its peak amplitude. 

Our results show that, as long as x is small, the dependence of the acceleration on the 
peak amplitude is in fact very similar to the one obtained by perturbation, that is, by using 
expressions ([S]) or ([7]) with y/C replaced by Ppeak- Let us first discuss the results for Pth = 0. 
Figure [1] shows the acceleration as function of peak amplitude of -P(C) for three different 
strengths of the plasma term and without the Raman term. As shown, the acceleration 
increases with the amplitude of the peak and a good agreement exists between the shooting 
results and the perturbation expression (jTj). Nevertheless, for x = 0.3 there is an observable 
difference between the two results that can be attributed to the deviation of the pulse profile 
from the sech shape considered in equation ([71). The absence of results for low Ppeak in the 
curves for x = 0.2 and x = 0.3 is due to the inability of our shooting to find a solution with 
a small right tail. In fact, the pulse profiles with small amplitude are located close to the 
zero of the Airy z axis, effect that is more pronounced as x grows. This means that, in those 
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FIG. 1. Dependence of acceleration parameter a on the pulse peak amplitude for Pth = and for 
three different x values. Points are shooting results and lines are for the perturbation expression. 





(a) (b) 

FIG. 2. Pulse profiles with peak amplitude close to (a) 0.3 and (b) 0.8 for two different x values. 

situations, the solution is no longer similar to a sech profile but it has long, and eventually 
oscillatory, tails. Figure |2] presents two sets of solutions for Ppcak ~ 0.3 and Ppeak ~ 0.8. The 
first set shows considerable differences at the right tail, with the pulse for x = 0.2 having 
a longer tail. In the second set, the shape differences are not so evident since for this peak 
amplitude, both solutions are already similar to each other and with the sech shape. 

Concerning the numerical results for Pth 7^ 0, the accelerations are again in good agree- 
ment with the perturbation expression if x is small. As figure 12] shows, when compared with 
the results for Pth = 0, the accelerations are lower for small peak amplitudes and larger 
otherwise. Furthermore, the acceleration decreases with increasing Pth for smaller peak am- 
plitudes and the inverse is true for larger peak amplitudes. Also represented in this figure is 
the acceleration resulting from the approximate expression ([8]) with \fC replaced by Ppeak, 
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FIG. 3. Dependence of acceleration parameter a on the pulse peak amplitude for different Pth 
values and fixed x = O-l- Points are shooting results and lines are for the perturbation expression 
([6]). The inset shows the shooting accelerations for Pth = 0.2 compared with a computed using ([8]) 
with ^/C replaced by -Ppeak- 
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FIG. 4. Pulse profiles with peak amplitude close to (a) 0.25 and (b) 0.5 for two different Pth values 
and fixed x = 0.1. 

which, as anticipated, is close to the shooting results when Ppeak ~ -Pth- The pulse shape 
differs from the sech also for some peak amplitudes close to Pth but this effect is smaller 
as Pth increases. Note that as Pth increases, the possible peak amplitudes that make the 
plasma term nonzero are also increasing, since the latter should be larger than the first. 
Figure m compares pulse profiles for Pth = and 0.2. 

As we introduce R different from zero, i.e., include the Raman term, the behavior of the 
acceleration with the peak amplitude is illustrated in Figure |5l For small peak amplitudes. 
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FIG. 5. Dependence of acceleration parameter a on the pulse peak amplitude for different R 
values and fixed x = 0.1 and Pth = 0. Points are shooting results and lines are for the perturbation 
expression. 

a increases with the peak amphtude whose behavior is characteristic of the plasma effect. 
However, as the peak amplitude increases further, the acceleration starts to decrease into 
the region of negative accelerations that are characteristic of the accelerating solitons of the 



NLS plus IRS [13|]. Note that, similarly to what happened when only the plasma term was 
present, the shooting was performed forward, but in the cases of negative a, the estimates 
of P and P' in the left tail were taken from Bi. Still for the case R ^ 0, the proximity of the 
pulses to the z = in the Airy axis for smaller peak amplitudes is lost as we introduce the 
Raman term. However, as the peak amplitudes increase to values for which the acceleration 
is negative and large in modulus, the pulse returns to the neighborhood of the Airy z = 
and starts to develop long tails but, in this case, to the left. 

Let us return to the case R = 0. Our results indicate that the pulse profiles develop 
long tails whenever the peak amplitude is small. Since these long tails can be associated 
with the oscillatory behavior of the Airy functions for negative z, we plotted in figure [6] the 
relative pulse amplitude at 2; = as a function of the peak amplitude for different values of 
X and Pth- As expected, this relative amplitude increases with x- On the other hand, when 
Pth = 0, this relative amplitude increases with the decrease of Ppeak but, for Pth 7^ the 
curves exhibit a maximum for a given value of Ppeak that is not large when compared to Pth- 
The existence of the long tails for small peak amplitude pulses is not readily understood 
since for those amplitudes the plasma term is smaller. Also, we know that in the presence of 
IRS the pulses develop long tails if the Raman term is large, which happens for large tr or 
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FIG. 6. Dependence of the relative amplitude at z = on the pulse peak amplitude for different 
X and Pth values. 
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(a) (b) 

FIG. 7. Effective nonlinear index for peak amplitudes close to (a) 0.5 and (b) 0.1 for fixed x = O-l- 
Comparison with the NLS case. 



short pulses (large peak amplitudes). In order to better understand the deviation from the 
sech shape in the presence of the plasma term, we plotted the effective nonlinear refractive 
index |gp — Tji(\q\'^)r — 0t ( 1 — e~°^'^-°° '^' "*• m ) '^ \ for several cases. Figure [7] compares 
two of those cases against the typical nonlinear refractive index of the NLS. We may observe 
that, although the magnitude of the effect of the plasma term is greater for -Fpeak = 0.5, the 
relative deviation from the NLS case is larger in the -Fpeak = 0.1 case. There, we may also 
see that the introduction of nonzero Pth decreases the plasma effect what was fully expected 
since this means that only part of the pulse creates the plasma. 
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FIG. 8. Evolution of pulse solutions, \q\ with (a) peak amplitude equal to 0.5 for x = 0-1 ^^^ 
Pth = 0.2 and (b) peak amplitude equal to 0.3 for x = 0.2 and Pth = 0. We have used a = 1. 

IV. DIRECT NUMERICAL INTEGRATION 



Direct numerical simulations of the full equation ([3]) were performed using pseudospectral 
codes in order to study the stability of the solutions described in the last section and to 
confirm accelerations and existence of tails. In general, the solutions found by the shooting 
procedure are stable and evolve along the predicted trajectory. Figure 8(a) is a contour 
graph showing the evolution of the pulse profile of Fig. 4(b) for Pth 



0.2. The trajectory 



is in full agreement with the predicted acceleration. Whenever long tails were found in the 
shooting procedure, they were confirmed in the propagating solution. In cases of very long 



tails, the solution is no longer stable but decays. Figures 8(b) and [9] show the evolution and 



output as obtained for the pulse profile of Fig. 2(a) corresponding to x = 0.2. The same 
kind of behavior was already observed with the Raman accelerating solutions [13] and is 
consistent with infinite energy of the Airy solutions Ai{z) and Bi{z). As discussed in section 
llllt pulse solutions in the positive Airy axis and far away from its zero, have exponential 
decay in both tails, similar to Ai to the left and to Bi to the right (the inverse happens for 
a < 0). The algebraically decaying oscillations of Bi(z) for negative z will only occur far 
way in the right tails (left tails for a < 0). However, if the solutions are in the positive Airy 
axis but close to the zero one of the tails will behave like a combination of Ai and Bi, it will 
exhibit the typical algebraically decay and it will shed radiation away during propagation. 



Figures 8(b) and [9] report these latter behavior. 

Finally, let us return to the physical variables and calculate the actual acceleration and 
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FIG. 9. Input and output for the simulation whose contour is on Fig. |8(b)[ 

frequency shift. The adimensional acceleration a does not correspond directly to an acceler- 
ation in physical units, nevertheless let us define the acceleration a^ as the second derivative 
of the temporal peak position in the group velocity reference frame relatively to the propa- 
gation distance z, namely 



d tpeak _ 1/32 1 d Tpcak _ |/32| Ct 



dz^' 



''0 



de 



tl 2 



2a3|/32 



-a(X,^th,^pcak) 



Note that the negative signal only implies that the pulse is traveling toward negative t but 
since t is measured in a reference frame that travels with the group velocity for uo, the pulse 
is gaining velocity whenever this acceleration is negative. This change in velocity is due to a 
deviation in frequency that is linear with the distance z, as expressed in the phase @, and 
given by 

d9 a^ 



Au 



Since equation (^ neglects the photoionization related losses that are small for pulses whose 
peak amplitude is comparable with the threshold, for x not too large, we may use expression 
([8]) for a and approximate the frequency approximated by 



Acu 



15 1/521 



^peak^ + 



V^peak 



3/2 



3A 



which gives the standard result for the IRS [10|, ll3[ and the effect of plasma growing with 
order ^2^^^- 
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V. CONCLUSIONS 

We have found the self-similar accelerating solutions of a generalized NLS that includes 
IRS and a term for plasma induced nonlinearity. This equation models the propagation of 
pulses in gas filled HC-PCFs where it has occurred photoionization of the gas. The solutions 
are very close to the NLS sech soliton as long as the strength of the plasma term is relatively 
low and the solution amplitude is relatively large. The accelerations and the blueshifting 
increase with the peak amplitude of the pulses. In case of pulse solutions whose peak 
intensity is close to the photoionization threshold, which are the ones for which the equation 
better models the physical effects, the frequency blueshift increases in the same order as 
the square of peak amplitude. However, also the same solutions, whose peak amplitudes 
are close to the threshold, may exhibit a profile that is considerably different from the sech, 
have long tails and decay along the propagation distance. 
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